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ABSTRACT 

We present a largely improved version of CRASH, a 3-D radiative transfer code that 
treats the effects of ionizing radiation propagating through a given inhomogeneous 
H/He cosmological density field, on the physical conditions of the gas. The code, 
based on a Monte Carlo technique, self-consistently calculates the time evolution of 
gas temperature and ionization fractions due to an arbitrary number of point/extended 
sources and/or diffuse background radiation with given spectra. In addition, the ef- 
fects of diffuse ionizing radiation following recombinations of ionized atoms have been 
included. After a complete description of the numerical scheme, to demonstrate the 
performances, accuracy, convergency and robustness of the code we present four dif- 
ferent test cases designed to investigate specific aspects of radiative transfer: (i) pure 
hydrogen isothermal Stromgren sphere ; (ii) realistic Stromgren spheres; (iii) multiple 
overlapping point sources, and (iv) shadowing of background radiation by an inter- 
vening optically thick layer. When possible, detailed quantitative comparison of the 
results against either analytical solutions or 1-D standard photoionization codes has 
been made showing a good level of agreement. For more complicated tests the code 
yields physically plausible results, which could be eventually checked only by compar- 
ison with other similar codes. Finally, we briefly discuss future possible developments 
and cosmological applications of the code. 

Key words: cosmology: radiative transfer - intergalactic medium - Lya forest - 
reionization 



1 INTRODUCTION 

Stimulated by the growing number of observations of the 
high-redshift universe {e.g. Fan et al. 2001, 2003; Hu et al. 
2002), an increasing attention has been dedicated to the the- 
oretical modeling of the reionization process, adopting both 
semi-analytic and numerical approaches {e.g. Gnedin & Os- 
triker 1997; Haiman & Loeb 1997; Valageas & Silk 1999; 
Ciardi et al. 2000, CFGJ; Miralda-Escude, Haehnelt & Rees 
2000; Chiu & Ostriker 2000; Gnedin 2000; Benson et al. 
2001; Razoumov et al. 2002; Ciardi, Stoehr & White 2003, 
CSW). An important refinement has been introduced by the 
proper treatment of a number of feedback effects ranging 
from the mechanical energy injection to the H2 photodisso- 
ciating radiation produced by massive stars (CFGJ). Only 
recently, though, the reionization process has been studied in 
a cosmological context including structure evolution, a real- 
istic galaxy population and a sufficiently accurate treatment 
of the Radiative Transfer (RT) of ionizing photons through 
the InterGalactic Medium (IGM), in a simulation volume 
large enough to have properties "representative" of the all 
universe (CSW). The next challenge for this kind of study is 



to develop accurate and fast radiative transfer schemes that 
can be then implemented in cosmological simulations. 

Although the reionization process is one of the most 
demanding and natural applications of the radiative trans- 
fer theory, numerous physical problems require a detailed 
understanding of the propagation of photons into differ- 
ent environments, ranging from intergalactic and interstel- 
lar medium to stellar or planetary atmospheres. The full 
solution of the seven dimension RT equation (three spatial 
coordinates, two angles, frequency and time) is still well be- 
yond our computational capabilities and, although in some 
specific cases it is possible to reduce its dimensionality, for 
more general problems no spatial symmetry can be invoked. 
Thus, an increasing effort has been devoted to the develop- 
ment of radiative transfer codes based on a variety of ap- 
proaches and approximations {e.g. Umemura, Nakamoto & 
Susa 1999; Razoumov & Scott 1999; Abel, Norman & Madau 
1999; Gnedin 2000; Ciardi et al. 2001, CFMR; Gnedin & 
Abel 2001; Gen 2002). In the following, we give a list of 
the available cosmological radiative transfer codes, together 
with their characteristics and applications. 

• Gnedin & Ostriker (1997): time-dependent code based 
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on the local optical depth approximation. It solves the evo- 
lution of H, He, H2 and temperature. It includes the RT 
of radiation from multiple point sources, a background flux 
and the diffuse component due to recombination. The code 
is fast, independent on the number of sources, but the ap- 
proximation breaks down for optical depths of the order of 
one and cannot account for shadowing effects. It has been 
applied to the study of the reionization process and has been 
coupled with hydrodynamical simulations in the implemen- 
tation of Gnedin (2000). 

• Umemura, Nakamoto & Susa (1999): time-independent 
code based on the method of the short characteristics. It cal- 
culates the equilibrium configuration of an isothermal H/He 
gas. It includes the RT of a background flux and the diffuse 
component due to recombination. As a backdraw the code 
can calculate only equilibrium configurations. It has been 
applied to the study of cosmological reionization. 

• Abel, Norman & Madau (1999): timc-dopcndcnt ray- 
tracing code. It solves the evolution of hydrogen and treats 
RT of the radiation emitted by a point source, together with 
the diffuse component. Although the ray-tracing approach 
is very accurate, its primary limitation is the high computa- 
tional cost required to follow a largo number of sources and 
to assure the correct coverage of the simulation volume. This 
code has been applied to the study of the evolution of the 
ionized region surrounding a mini-quasar in a cosmological 
density field. 

• Razoumov & Scott (1999): time-dependent ray-tracing 
code, based on the solution of the 5D (three spatial coor- 
dinates and two angles) advection equation. It includes the 
H, He and H2 chemistry and the temperature evolution. It 
treats the RT of radiation from multiple point sources, a 
background fiux and the diffuse component. The code suf- 
fers from the same problem as the previous one. It has been 
applied to the study of the reionization process. 

• Ciardi et al. (2001): CRASH 1.0 time-dependent code, 
based on a Monte Carlo approach. It includes the evolution 
of an isothermal H gas and it treats the RT of the radiation 
emitted by point sources, together with the diffuse compo- 
nent. The code has been applied to the study of the evolution 
of the ionized region surrounding point sources embedded in 
a cosmological density field. 

• Gnedin & Abel (2001): time-dependent code, based on 
the exphcit moment solver OTVET (Optically Thin Vari- 
able Eddington Tensor formalism). It includes the H and 
He chemistry and the temperature evolution. It treats the 
radiation produced by point sources, a background flux and 
the diffuse component. Its main problem is related to the 
correct evaluation of the Eddington tensors. In fact, while 
the assumption of optically thin Eddington tensor is exact 
in the limit of a single point source, it may break down in 
complex situations with very inhomogeneous source func- 
tions and opacity fields. It has been coupled with a hydro- 
dynamical simulation to study cosmological reionization. 

• Sokasian, Abel & Hernquist (2001): time- dependent 
ray-tracing code. It includes the H and He chemistry and 
the temperature evolution. It treats the radiation produced 
by point sources, a background field and the diffuse compo- 
nent. As a first application, it has been used to the compute 
Hell reionization by quasars. 

• Cen (2002): time-dependent ray-tracing code, based 
on the fast Fourier transform method for hydrogen ioniza^ 



tion. It treats the radiation produced by point sources, a 
background field and the diffuse component. Unlike conven- 
tional ray-tracing schemes where angular discretization is 
performed on the source, here it is done at the receiving 
site, guaranteeing good space coverage regardless of the fi- 
nesse of the angular discretization. Moreover, its speed is 
independent on the number of sources. The code has been 
applied to test cases of single and multiple sources embedded 
in a constant or inhomogeneous density field. 

• Razoumov et al. (2002): time-dependent ray-tracing 
code. It includes the H, He and H2 chemistry and the tem- 
perature evolution. It treats the radiation produced by point 
sources, a background field and the diffuse component. To 
alleviate the speed problem related to the ray-tracing ap- 
proach, the authors have implemented the method devel- 
oped by Abel & Wandelt (2002), introducing trees of rays 
segments which recursively split into sub-segment as one 
moves away from the source. The code has been applied to 
the study of the reionization process. 

Monte Carlo (MC) methods, on which we base our ap- 
proach, have been widely used in several (astro-)physical ar- 
eas to tackle RT problems (for a reference book see Cashwell 
& Everett 1959) and they have been shown to result in fast 
and accurate schemes. Here we build up on previous work 
(Bianchi, Ferrara & Giovanardi 1996; Ferrara et al. 1996; 
Ferrara et al. 1999; Bianchi et al. 2000) to develop a new 
version of the code CRASH [Cosmological Radiative transfer 
Scheme for Hydrodynamics) , which improves in many ways 
the earlier version presented in CFMR. The main aim of 
this paper is to give a full description of the code CRASH in 
its present form, perform a number of stringent tests that 
demonstrate its ability to treat RT problems arising in a 
wide range of astrophysical context, and evaluate its per- 
formances. In forthcoming work, we will apply the code to 
specific cosmological problems. 

The plan of the paper is as follows. In Sec. 2 we de- 
scribe the details of the numerical implementation concern- 
ing the discretization of the radiation field and its (time- 
dependent) interaction with the surrounding matter; we pay 
special attention to issues as the dependence of spatial and 
time resolution on the numerical parameters. In Sec. 3, we 
present results from a set of specifically designed test-runs: 
Stromgren spheres, two point sources and an optically thick 
slab case illuminated by a background field. Finally, in Sec. 
4, we draw our conclusions and discuss some of the possible 
future applications and developments of CRASH . 



2 DESCRIPTION OF THE CODE 

CRASH , is a 3-D numerical code primarily developed to study 
time-dependent radiative transfer problems in cosmology; 

however, the versatility of its scheme allows easy extensions 
to additional astrophysical applications. Such a scheme is 
largely based on MC techniques, which are used to sam- 
ple appropriate Probability Distribution Functions (PDFs), 
as briefly explained in Appendix A. The current version of 
the code (CRASH 2.0) is a development of a previous imple- 
mentation of the same scheme, discussed in CFMR. The 
main differences with that first attempt are [i] the inclusion 
of helium in all its ionization states, [ii] the self-consistent 



CRASH: a Radiative Transfer Scheme 3 



calculation of gas temperature, and [iii] an efficient algo- 
rithm to treat multiple radiation sources and/or a diffuse 
background field. In addition, more refined photon propa- 
gation schemes (see below) have been included. As today, 
CRASH works as a stand-alone routine that can be applied to 
arbitrary and pre-computed hydrodynamical density fields. 
A future development will consist of a full coupling with 
available cosmological hydrodynamical codes, so that the 
interaction of radiative transfer and gas dynamics can be 
treated sclf-consistcntly. Wc now turn to the description of 
the implementation of the physical processes and the nu- 
merical algorithms used by the code. 



2.1 Initial Conditions 

As a first step, the initial couditious for the physical quan- 
tities (gas density, temperature, H and He ionization frac- 
tions) must be assigned on the iV^ nodes of the adopted 
3-D cartesian grid. The density field, the temperature and 
the ionization fractions can be either the output of a hydro- 
dynamical simulation or arbitrarily constructed. While the 
density field remains constant during the simulation, i.e. no 
back-reaction of gas dynamics to the effects of the radia- 
tive transfer is considered, the temperature and ionization 
fractions are instead updated with time. Both open and pe- 
riodic boundary conditions can be adopted. In the latter 
case, the ionizing radiation that exits with a given direction 
from a cell {xf,,yb,Zb) on the boundary surface of the sim- 
ulation box, is rc-cmittcd with the same direction from the 
cell {—Xb,—yb,—Zb); the origin of the cartesian coordinate 
system is located at the box center. 

Finally, the properties of the radiation field must be 
specified. The quantities required arc the number, location 
and emission properties (e.g. the intensity of the emitted rsr 
diation and its spectral energy distribution) of point sources 
in the simulated volume, and the intensity and spectral en- 
ergy distribution of the background radiation, if present. As 
CRASH has been specifically designed to follow the propaga- 
tion of ionizing photons, the above spectral parameters are 
required only for energies > 1 Ryd. The required proper- 
ties of the radiation fields are described in detail in the next 
Section. 



2.2 Ionizing Radiation Field 

The MC approach to radiative transfer requires that the ra- 
diation field is discretized into photon packets; the proper- 
ties of the radiation field and the relevant radiation-matter 
interaction processes are then statistically treated by ran- 
domly sampling appropriate PDFs. In this approach the 
ionizing radiation emitted by either a specified number of 
point sources located arbitrarily in the box, or a background 
radiation, is easily reproduced. Also, difi'usc radiation from 
recombinations in the ionized gas can be treated in a similar 
manner without further assumptions. 



2.2.1 Point Sources 

Let us first consider the case of a single point source with 
a time-dependent bolometric luminosity Ls{t). If ts is the 



(physical) simulation time, the total energy emitted by the 
source is: 

E,= f Ls{t)dt. (1) 
Jo 

Such energy is distributed into A'p photon packets, emit- 
ted at the source location at regularly spaced time inter- 
vals dt — ts/Np. Thus, the time resolution of a given run 
is fixed by Np. To each emitted photon packet a frequency, 
u, and a propagation direction, {0,4>), are assigned by sam- 
pling the source spectrum, Su, and the angular PDF, using 
the MC method described in Appendix A. This procedure 
allows to treat any given spectrum and possible emission 
anisotropics. The j-th monochromatic packet of frequency 

1/ is then emitted at time tj = jdt (j = 1, ,Np), and it 

contains N^,,j = AEj/hv photons, where AEj is the energy 
of the j-th packet given by: 

AEj = I ' Ls{t)dt. (2) 

Once the quantities v, {8, (f>) and N^j are assigned to an 
emitted photon packet, we follow its propagation from the 
source location into the given density field as described in 
Sec. 2.3. 



2.2.2 Background Radiation 

The presence of an external ionizing background radiation 
is treated with an algorithm analogous to the one used for 
point sources, with the only difference that the emitting lo- 
cation is one of the 6 x cells on the faces of the simulation 
box, rather than the cell corresponding to the source loca- 
tion. If the background radiation is uniform and isotropic, 
an integer number is randomly extracted in [1, 2, 3, 4, 5, 6] 
to select a face of the box and two random integer numbers 
in [1, Nc] to select the cell position on the chosen face. This 
procedure can be easily generalized to account for possible 
asymmetries in the radiation field. 



2.2.3 Diffuse Radiation 

Diffuse ionizing radiation in the IGM can be produced by 
recombinations and by free-free (f-f) transition processes. 

Nevertheless the f-f cmissivity at frequencies higher than 
1 Ryd is negligible with respect to the emissivity due to 
recombinations and we can safely exclude this process. The 
implementation includes as sources of diffuse radiation the 
recombination processes occuring in a H/He gas: 

• H+ -I- e" — > H° + hi^, 

• He+ + e- — > He" + hiy, 

• He++ + e" — > Hc+ + hi^. 

The MC scheme allows a straightforward and self-consistent 
treatment of the diffuse radiation produced by H/He recom- 
binations in the ionized gas. In analogy with the direct one, 
the diffuse radiation field is discretized into photon packets 
whose propagation is followed using the scheme discussed in 
the following Section. For this approach, it is necessary to 
assign the angular and frequency distribution of the emitted 
diffuse radiation as well as its intesity. 
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As the emission process for diffuse radiation is isotropic, 
the packet propagation direction is randomly selected. The 
spectral distribution of the diffuse radiation depends on the 
emitting species. If LTE is assumed, the emissivity rjii/) as- 
sociated to an arbitrary recombining atom, results in the 
following frequency dependence (Mihalas 1978, Osterbrok 
1989): 



> i'th,n, (3) 



nn{v) oc aH(t/)i^^e-<'"'-'"'*''.«''"'«^ 

where aniy) and huth.n are the photoionization cross- 
section and the ionization potential of the recombined atom 
H; kB is the Boltzmann constant and T is the kinetic tem- 
perature of the recombining electron. In the calculation we 
include only the effect of the ionizing diffuse radiation. 

The intensity of the diffuse radiation is evaluated as 
follows. For each of the three species we allocate a 3-D ar- 
ray, Nrec{x,y, z), composed of A*'^ elements, i.e. one for each 
cell of the 3-D cartesian grid. Each element of Nrcc(x,y, z) 
keeps track of the number of recombination events occurred 
in the corresponding cell. For simplicity, we now consider 
a fixed cell and the recombination process of species I 
[I G {H"^, He"*", He"*"^}) to describe how Nrec{x, y, z) is com- 
puted. If a packet crosses the fixed cell at a time tc, the 
number of recombinations, AA^rec, occurred since the last 
time, tc-i, a packet crossed the same cell is approximately 
given by: 

^3 



ANr. 



a{Tc-i)ne,c-ini,c-iAt{Ax) 



(4) 



where Tc-i,ne,c-i and n/,c-i are the temperature, electron 
number density and number density of species / evaluated at 
time tc-i; a{Tc-i) is the total recombination coefficient of 
species I, At — tc — tc-i is an integer multiple of dt, and Aa; 
the linear dimension of a cell. We then update the number 
of recombinations in the cell, Nrec.c = Nrec,c-i + AA^rec. If 
Nrec,c satisfies the following condition: 



a packet containing A'^^ 



A-,, 



(5) 

photons is emitted from 



the cell and Nrec,c is set to zero; we assume that all Nrec,c 
recombinations occur to the same atomic level, so that 
the emitted packet is monochromatic. In the above equa- 
tion, fr € [0, 1] is an adjustable numerical parameter and 
Na is the total number of recombining atoms in the cell, 
i.e. naiAx)^ or nHc(Aa;)^, where nn ~ nno + nH+ and 
riHe = riHeO + n-^^+ + nHe++ ■ In principle low values of fr 
make the simulated recombination process smoother at the 
expense of a higher computational cost; experimentally we 
have found that /,. =0.1 represents a good compromise be- 
tween accuracy and computational cost. 

At each recombination event, e.g. each time the condi- 
tion |5] is satisfied, we stochastically derive the probability 
for recombinations at ground level to occur from the ratio 
Qi(r)/a(r), with Qi being the recombination coefficient to 
the ground level (the expressions for all the rates adopted in 
CRASH are given in Appendix B). If recombination occurs to 
the first level we sample the spectrum in eq.|21 (appropriately 
normalized) for all the three reemission processes. Recom- 
binations to levels other than the first of H'' are neglected 
as they produce no ionizing photons. De-excitations follow- 
ing recombinations to higher levels of He" or He"*" produce 
instead ionizing photons for H" or He°, respectively. The 
de-excitation can occur through different channels; however 



for this cases the frequency assigned to the emitted packet 
is the one corresponding to the relevant transition from the 
second to the ground level: for He" we assume the mean 
value between hv = 19.8 eV and hv = 21.2 eV (2s^ls and 
2p^ls); for He^ hu — 40.7eV (2p^ls). This approximation 
has negligible effects on the results of a typical cosmological 
simulation and it can be easily modified if a more detailed 
treatment of chemistry is necessary. 



2.2.4 Multiple Sources 

As already pointed out, CRASH allows to run simulations with 
more than one source. Let us consider the case of A^s sources 
(one of which could be a background radiation field) emit- 
ting simultaneously ionizing radiation inside the box. Each 
of the Ns sources is treated as described in the previous 
Sections, assuming the same value of Np for all of them. In 
this way it is possible to reproduce the time evolution of the 
total ionizing radiation field in the computational volume, 
by emitting the j-th packet {j £ {l,..,A'p}) from all the 
Ns sources at the same time-step tj = jdt, looping on the 
sources at each j value. 

This simple method has been chosen as it is computa- 
tionally very cheap. We have already noted that the time 
resolution and the accuracy (and the convergency, see later 
on) of the simulation depends on the total number of pack- 
ets emitted by a single source, Np. In general, to achieve 
the same accuracy level when multiple sources are present, 
one would need to multiply such number by the number of 
sources in the box, thus becoming NpNa, with a significant 
increase of the computational time. However, in practice this 
is not necessary because, unless either the ionized regions 
created by each source are completely separated and/or the 
ambient gas is very optically thick, almost inevitably a typi- 
cal cell in the computational domain will "see" a large num- 
ber of sources, hence the effective number of photon packets 
required is ^ NpNa. 



2.3 Photon Packets Propagation 

The MC approach to the RT problem is particularly conve- 
nient. In fact, as the radiation is modeled in its particle-like 
nature, it is not necessary to solve the high-dimensional cos- 
mological RT equation for the specific intensity of the radia- 
tion, which usually requires several analytical and numerical 
approximations (Abel et al. 1999; Razoumov & Scott 1999; 
Nakamoto et al. 2001; Razoumov et al. 2002). 

To describe the propagation direction of the emitted 
packets we adopt spherical coordinates, (r, 9, 0), with origin 
at the emission cell, {xe,ye, Ze). Thus, a packet will propa- 
gate along the direction identified by: 



X = Xc + r sin 6 cos < 
y = ye + r sin sin <; 
z = Ze + r cos 9 



(6) 



Here, {9, (f>) are determined via the MC method applied to 
the appropriate emission angular PDF of the source (see 
Sec. 2.2.1). The cells crossed by the packet can be ordered 
by the index / = Q, ...,lmax (where Imax is the maximum 
number of cells that are crossed by the packet), which is 
related to the cells cartesian coordinates by the following 
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expression: 

x{l) = Xe + int{l ao) 

y(l) = tJe + int(l bo) 
z{l) — Ze + int{l Co). 



(7) 



Here ao , 60 and co are the director cosines of the propagation 
direction with respect to the cartesian grid: ao = sin 9 cos 4>, 
bo = sin ^ sine/), cq = cos^. 

Given the above definitions, let us now focus on the 
propagation of the generic monochromatic packet initially 
composed of A''-y photons of frequency i/. As already men- 
tioned, photoelectric absorption of continuum ionizing pho- 
tons is the only opacity contribution included. Thus, the ab- 
sorption probability for a single photon travelling through 
an optical depth r is given by: 



P(r) 



(8) 



As the packet propagates through the above cell sequence, 
it will deposit in the l-th cell a fraction of its initial photon 
content oc P{t^), where r' is the total optical depth crossed 
from the emission site to the cell itself. For each cell crossed, 
we compute the contribution, A'r, of the actual cell to the 
total optical depth, as follows: 

A't = AV^o + A'r^go + AV, 



[a^o(i/)n^o + 
cr^<jo(i/)n^^o +CTHe+(i')'"4e+]/W^a;, 



J/e+ 



(9) 

where ga is the pliotoionization cross-section for absorber 
A £ He", i/e+}, n/4 its numerical density in l-th cell 

and f{l)Ax the path length through the Z-th cell of linear 
size Aa;. Depending on the trajectory of the packets in the 
cell, f{l) will have a value in the range [0,V3]. To limit 
the computational cost, wc do not include in the code a ray- 
casting routine to determine the path length in each cell. We 
assume, instead, the fixed value /(Z) = 0.56, corresponding 
to the median value of the PDF for the lengths of randomly 
oriented paths entering from the faces of a cubic box of unit 
size (see CFMR). 

The number of photons absorbed in the Z-th cell, N^, is 
determined as follows. When a packet reaches the Z-th cell, 
its photon content is: 



■A'. 



A;-^e-^ 



" < A, 



(10) 



it will deposit in the same cell a number of photons equal 
to: 

Ai = A^(l - e""^'^). (11) 

In some cases can exceed the total number of atoms 
that can be ionized in the Z-th cell, A^^^; when this occurs 
wo sot = Nlg,-^ and retain the excess photons in the 
traveling packet. Propagation of each packet is followed until 
it exits from the box (open boundary condition case) or 
until extinction occurs when A^ < 10~^At,. This procedure 
guarantees energy conservation to an equivalent accuracy 
of 10"''. We typically adopt p € [4,9] depending on the 
required accuracy. 



2.4 Updating Physical Quantities 

In this Section we describe the scheme adopted to calcu- 
late the time evolution of the ionization fractions (a;^+ = 



Uu+friH, = nu^+friHe, Xh^++ = nu^++ /uHe) and 

of the gas temperature, T. The system of coupled equations 
that describes the above evolution is the following: 



riH 



riHe- 



dxjj+_ 
dt 

d xue+ 
dt 



dx 



He++ 



dt 



7^o(r)n^one - aH+{T)nH+ne 

7HeoiT)nH^one - jHe+iThne+ne - 

+ ^He°n„^o = Tu^+ , (12) 
'yHe+(T)nH^+ne - aH^++{T)nne++ne 



dT 

dt 



3fcBn I dt 



A{T,xi) 



The last equation is the energy conservation equation, where 
n = riH + riHe + ne is the number of free particles per unit 
volume in the gas; Ti. and A arc the heating and cooling 
functions which account for the energy gained and lost from 
the gas in the unit volume per unit time, respectively (see 
Appendix B for the detailed A expression adopted). In the 
three equations for the ionization fractions we indicate with 
Qj (ta) the recombination (coUisional ionization) coefficient 
and I e {H+, He+, He++} {A £ {H", He", He+}); Fa is the 
time- dependent photo-ionization rate. 

The numerical approach requires to discretize the dif- 
ferential equations as follows: 



Xh+ (t + At) = Xn+ it) -h Th+ {t)At/nn , 

XH,+ {t + At) = XHe+it)+THe+{t)^t/nHe, 



XH^++{t + At) 

T{t + At) 



XHe++{t) +lHe++{'t)^t/nHe, 

At[n{T,xi)-A{T,xi)]}, 



(13) 



where An = n{t + At) — n{t). The physical quantities in a 
cell arc updated by solving the above system each time the 
cell is crossed by a packet. Thus, the integration time step 
At {i.e. the time interval between two subsequent passages 
of a packet in a given cell) is not constant, due to the statis- 
tical description of the emission and propagation processes. 
The integration time step is calculated as follows. We in- 
troduce a 3-D array whose A^ elements correspond to the 
index of the last packet that has crossed a cell (Sec. 2.2.1). 
Thus, when the j-th packet crosses the cell with coordinates 
{x,y,z) wc update the physical quantities in the cell, using 
an integration time step equal to: 

At^[j^j'{x,y,z)]dt, (14) 

where the array element j'{x,y,z) gives the index of the 
packet last transited through {x,y,z), and we update the 
array element to the actual index j. 

As radiation field is discretized in photon packets, it is 
not straightforward to recover continuous quantities that di- 
rectly depend on its intensity, as photo-ionization and photo- 
heating rates. Their effects are evaluated by means of dis- 
cretized contributions as a function of the number of pho- 
tons deposited in the cell, N'a = A^_^o + Ajj^^^o + Aj^ ^^+; 
these are distributed among the three absorbing species pro- 
portionally to their absorption probability (1 — '^^). 
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The corresponding ionization fractions and temperature in- 
creases are given by: 



riH 
N\ ( 1 



A, go 

TV' 



n'-jjA^x \ 1 — e 



-At' 



Aa; 



He+ = 



1 Af : 



N' 



(15) 



Aa; 



'He 



He 



1 4- 
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Recombinations, collisional ionizations and cooling are in- 
stead treated as continuous processes. This requires At 
to be much smaller than the characteristic timescales 
of these processes for all species I and A, i.e. At <^ 
train = min{trsc,i ,tcoii.A,tcooi} ■ If this coudition is not 
fulfilled, the integration is split into steps, with Us = 
int [At/{fatmin)\, fs IS a fudgc factor, usually taken equal to 
fs = 50— 100, in order to minimize the discretization errors. 
The expressions used for recombination, collisional ioniza- 
tion and cooling rates are given explicitly in Appendix B. 



2.5 Numerical Resolution 

The input parameters that determine the spatial/time res- 
olution and the computational cost of a simulation are the 
following: 

• A'j,: number of photon packets emitted per source. 

• N^: number of grid cells. 

• A'^s: number of sources in the simulation. 

Although physical processes in our scheme are treated sta- 
tistically, it is nevertheless possible to derive semi-empirical 
dependences among these parameters and the performances 
of CRASH in terms of resolution and speed. Once the initial 
conditions are fixed, Np determines the mean number of pho- 
tons contained in a single packet, A'^, (actually the photon 
number content in each packet depends on its energy and 
frequency), and the time interval between two subsequent 
packet emissions, dt. The number of grid cells, A^^) deter- 
mines the spatial resolution and the mean number of atoms 
in each cell, Na. 

Let us consider the general case of a simulation with A'^s 
sources and open boundary conditions. The total number of 
mitted packets is Np°^ = NsNp. To a first order of magnitude 
a packet will pierce a mean number of cells equal to fdNc, 
where fd is a parameter ranging in [A^c~^, 1]: /d ~ 1 for an 



optically thin medium and fd ~ Nc~'^ for a simulation box 
with optically thick cells. During a simulation a cell will be 
crossed and updated an average number of times equal to 



N'i 



(16) 



the accuracy of a given simulation is mainly determined by 
the magnitude of this parameter. In order to correctly re- 
produce the physics of photoionization, in fact, the value 
of Ncr must be large enough to appropriately sample the 
frequency and the time distribution of the ionizing radia- 
tion at each cell location. This implies that it is necessary 
to impose a lower limit on N^r. We find that most common 
{e.g. power-law, black-body) spectra can be well MC sam- 
pled with roughly 10^ — 10* extractions {i. e. packets crossing 
each cell), depending on the shape of the spectra and on the 
required accuracy. Hence we require that Ncr exceeds the 
previous estimate, to achieve a good spectral sampling in 
each cell. A further condition has to be imposed on Ncr, 
which determines also the mean value of the integration 
time-step, {At) = ts/Ncr (see above). As already pointed 
out, At must be much smaller than the characteristic time 
scales of the physical processes treated as continuous. This 
imposes an upper limit to At, or, in turn, a lower limit to 
Ncr set by: 
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NsNp 
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(17) 



We initialize tmin = min{trcc,i, tcoii.A,tcooi} before the run, 
with the characteristic times calculated using the mean val- 
ues of the initial temperature and ionization fractions fields; 
although this quantities difi'er from cell to cell and will in- 
evitably change due to the evolution of the physical prop- 
erties of the gas during the simulation. However, if the con- 
dition eq. 1171 is not fulfilled at some computational stage, 
we resort again to splitting as described in Sec. 2.4. Before 
running a simulation, the numerical parameters have to be 
chosen such to satisfy the above conditions; then the accu- 
racy of the calculation increases with Ncr. 

Two other combination of the input parameters infiu- 
ence the accuracy of the calculation. The first one is the ratio 
N^/Na in a single cell which determines the mean number 
of packets necessary to completely ionize a single cell; the 
lower its value, the better the adopted scheme is able to 
reproduce the continuity of the photoionization (photoheat- 
ing) process. 

The second combination is the one giving the ratio be- 
tween the photoionization rate and that of continuous pro- 
cesses {e.g. recombinations) averaged in the simulation box. 
Taking into account the geometrical dilution of the ioniz- 
ing radiation emitted by point sources, a rough estimate for 
such quantity is < N-y > {Lbox)~ tmin, where < Nj > is the 
mean ionizing photon rate. As discussed in the previous Sec- 
tion, we adopt a different method to describe the effects of 
these two classes of physical processes. The scheme adopted 
to account for photoionization is very precise, as it guar- 
antees energy conservation to an adjustable accuracy 10"''. 
Recombinations, collisional ionizations and several other ra- 
diative processes which contribute to gas cooling, are instead 
reproduced as continuous by introducing the corresponding 
rate coeflicients. For this reason errors on these quantities 
are introduced by the discretization of the equations. Hence 
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Figure 1. Time evolution of the residuals of the Stromgren sphere 
radius, 7?,^, for a run with A^p = 10^ and y=6,7,8, calculated with 
respect to = lO'' (see text for details). 

the level of accuracy can be controlled very well if photoion- 
ization is dominant with respect to continuous processes by 
simply increasing Nc.r\ if < ^Vt- > {Lijox)^^tmin becomes too 
low, the accuracy can instead be degraded: however, the evo- 
lution of this quantity is monitored in order to prevent such 
degradation to occur during the simulation. 
The computational cost of a given simulation is proportional 
to N^Ncr = fdNcNpNs, which gives a rough estimate of the 
number of times the system eg. I13l is solved. 



3 TESTS 

In this Section, we present a set of tests to illustrate the per- 
formances of CRASH in a wide range of astrophysical and cos- 
mological contexts. Radiative transfer problems in three di- 
mensions do not have, in general, analytical solutions except 
for some simple cases; for this reason testing quantitatively 
a radiative transfer code is a challenging purpose. As a first 
test we consider the case of a monochromatic point source 
embedded in a pure- hydrogen homogeneous gas, with the 
ionized component at constant temperature. This idealized 
case of an isothermal Stromgren sphere, has the advantage 
of being directly comparable to an analytical solution of the 
problem {e.g. Spitzer 1978). Hence, it is possible to verify 
the accuracy of the code in reproducing the time evolution of 
the simulated system and its performances in a large range 
of gas densities. A second test is designed to check the reli- 
ability of the algorithm in calculating the temperature and 
ionization structure of a gas of primordial H/He composi- 
tion. As no analytical solution is available for this more re- 
alistic case, as a check, we compare our results with the ones 
of the 1-D RT code CLOUDY94\ available on the web. The 
above test is extended to the case of two point sources with 
largely different luminosities. This problem turns out to be 
a difficult benchmark for RT codes (see Gnedin & Abel 2001 

^ http:/ /nimbus. pa. uky.edu/cloudy/ 
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Figure 2. Neutral hydrogen fraction, function of the 

distance from the source in cell units. The points represent a 
subsample of grid points for runs with: Np = 10®, 10^, 10*, 10^ 
(triangles, crosses, asterisks and squares). 

for a discussion). Finally, we study the case of a dense, very 
optically thick slab illuminated by a power-law background 
radiation: our main aim is to show how CRASH performs when 
confronted with the problem of the diffuse radiation field 
produced by recombining gas (shadowing). 
To test even more complex problems (i.e. time evolution, 
inhomogeneous density field) it is instead necessary to com- 
pare results obtained with different RT codes; such an at- 
tempt (dubbed TSU^) is currently under development^. 

3.1 Pure Hydrogen, Isothermal Stromgren Sphere 

As mentioned, we first test the code against the analyti- 
cal solution for the time evolution of the radius of the Hn 
region produced by a monochromatic source with constant 
ionizing rate and expanding in a homogeneous medium. The 
comparison requires an helium abundance equal to zero and 
a constant temperature in the Hn region; we fixed it equal 
to T = 10" K. The (homo geneous) density field in the sim- 
ulation box is initialized with a hydrogen number density 
of uh ~ 1 cm~^; the cartesian grid, which has a linear size 
Lbox = 70 pc, is composed of A'^^ = 128"^ cells. This density 
value has been chosen in order for the diffuse (recombina- 
tion) radiation effects to be evident, so to test the corre- 
sponding implementation. The monochromatic source emits 
photons with energy equal to 13.6 eV, is located at the cen- 
ter of the box and has an ionizing photon rate Nj = 10"*s~^. 
The ionized hydrogen fraction is initialized at coUisional 
equilibrium at T = 10* K, i.e. = 1.2 x 10"^. To ensure 
that equilibrium configuration is achieved the simulation is 
carried out for a physical time ts — 6 x lO"" yrs, roughly 
equal to five recombination times. 

The Hii region equivalent radius, R„ , can be derived as 
Rn = {3Vn/4ny^^ , where Vn is the volume of the ionized re- 
gion. The contribution of each cell {x, y, x) to Vn is assumed 

^ |http://www.arcetri.astro.it/science/cosmology/Tsu3/tsu3.html| 
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Figure 3. Left: time evolution of tlie numerical Hu region radius Rn (dashed lines) compared with the analytical solution Ra (solid). 
The curves refer to eight different values of nu = 10~' cm~^ and are labeled with the index, q, of the corresponding density. Bight: 
percentage errors on the numerical radius with respect to the analytical one for the same runs (and labels) as in the left panel. 



to be a fraction of the cell volume equal to its ionization 
fraction x-j^+(x,y, z): 



Xjj+{Ax) 



(18) 



As a single cell could be brought to ionization equilibrium in 
several time steps, this approach guarantees that the volume 
of the Hii region is not overestimated as the contribution 
of the partially ionized cells at the edge of the I-front is 
included. 

To check for convergency, we run four simulations with 
increasing Np. Fig. shows the residual, TV , of a given run 
with respect to the highest resolution run (A'p = 10^), de- 
fined as; 



(19) 



where y = 6,7,8 refers to the run with A'j, = 10^. Good con- 
vergency is achieved already for Np = 10^, when the residual 
is <7%. 

In Fig. 121 the equilibrium neutral hydrogen fraction, x-f^o = 
rij^o/nn, profile is shown for the four runs. In this case, 
a higher number of photon packets (A^p 10*) is needed 
to reach a good level of convergency. This result suggests 
that applications which require the determination of the vol- 
ume ionization structure (as cosmic reionization) are less 
demanding than problems that require the correct value 
of neutral hydrogen fraction {i.e. Lya forest simulations). 
Thus, depending on the application and the desired accu- 
racy, the required number of photon packets might vary; 
anyway CRASH converges for a reasonable number of photon 
packets. As A'j, = 10^ already gives good convergency of the 
calculated numerical radius to be compared with the ana- 
lytical solution, in the following applications we will adopt 
this value, unless otherwise stated. 

The time evolution of R„ can now be compared with 
the classical analytical solution (Spitzer, 1978): 



where Rs — (3A'.^/47rn|f as) is the Stromgren radius 
and as is the hydrogen recombination coefficient to lev- 
els higher than the first. The comparison of our numerical 
results for i?„ against the analytical solution, has been per- 
formed for eight different values of the density in the range 
n — 10~® — 10 cm~^. In each run A'j,, Nc and A'^ are kept 
constant, as specified above. The source has been located 
at a corner of the simulation box (of linear size equal to 
6-Rg/5), to maximize spatial resolution. With the physical 
simulation time set to ts = 51^^^, the equilibrium configura- 
tion is reached within 0.25% of the total computational time, 
according to eq. I20I I. As Ncr is fixed for all the runs, the 
accuracy itself is the same, as far as the photon packet distri- 
bution and the discretization errors are concerned. Similarly, 
the ratio N^/Na, which determines the accuracy to which 
the continuity of the photoionization process is reproduced, 
is the same for all the runs. In fact. 



N^JNa = (A'^t.A',r')/(nHAa:') 

-1 / r / -2/3n31 i 

oc riu I nHXJijj ) = const. 



(21) 



Ra(t) = Rs(l 



(20) 



Finally, the ratio between the ionizing photon emission 
rate and the rate of physical processes reproduced as 
continuous (recombinations, collisional ionizations etc.) in 
a given cell increases with density: N-y{Li,ox)~^tmin oc 
Nj{n'^^^)~^njj^ oc n^^. Hence, if Nj remains constant, 
the accuracy of the results decreases towards lower densities 
for the reasons explained in Section 2.5. 

This conclusion is supported by the numerical experi- 
ments shown in Fig. 13 showing the evolution of the numer- 
ical radius, Rn (dotted line), compared with the analytical 
solution, Ra (solid line), for eight values of the density con- 
sidered. The agreement is remarkably good. For densities 
> 10~*cm~^ errors remain through the entire simulation 
^ 3%; for lower densities they increase, according to expec- 
tations following the previous discussion, up to ~ 10% for 
71 jf — lO^'' — 10"'^ cm"''. We point out that these runs have 
the minimum resolution required to achieve convergency and 
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that accuracy can be easily improved using a larger value for 

These tests, embracing 8 dex in density, demonstrate 
that CRASH is quite robust and accurate over a wide range 
of physical conditions. 

3.2 Realistic Stromgren Spheres 

Next, we test the CRASH scheme for the calculation of tem- 
perature and helium ionization. As no analytical solution is 
available, we compare the results with those obtained with 
the public code CLOUDY94, a 1-D radiative transfer code 
which calculates equilibrium configurations of the various 
species. The physics included in the two codes is slightly 
different. CLOUDY94 treats atoms and ions as multilevel 
systems, while we consider all species as one-level systems; 
furthermore CLOUDY94 includes some physical processes, 
like heat conduction and gas pressure effects, which are 
not included in CRASH . We consider a point source, emit- 
ting as a black body at T = 60000 K with a luminosity 
L = 10^^ erg s~^, located at the center of the simulation box 
and embedded in a homogeneous density field with n = 1 
cm~'^ composed of hydrogen (90 % by number) and helium 
(10 %). The gas is initially completely neutral and at a tem- 
perature T = 10^ K in the entire simulation box, whose 
linear dimension has been fixed to Ltox = 128 pc. As an 
accurate determination of the neutral fractions and tem- 
perature inside the ionized regions requires A'p 10* (see 
previous Section), we adopt A'p — 10*. Moreover, spheri- 
cal symmetry is required to degenerate our 3-D solution to 
the 1-D geometry adopted by CLOUDY94. The compari- 
son is performed at a time ts = 6 x 10^ yr « 5trec (where 
t^ec is the characteristic time scale for hydrogen recombi- 
nation to levels other than the first), when the equilibrium 
configuration has been reached. In Fig.|l|the comparison be- 
tween the results of the CLOUDY94 (sohd lines) and CRASH 
(points) is shown. The value of the different physical quanti- 
ties is plotted as a function of the distance from the source, 
expressed in cell units. Ax = 1 pc. The points represent 
spherically averaged CRASH outputs. Despite the differences 
between the two codes, a remarkably good agreement is ob- 
tained. In particular, the temperature profiles agrees within 
10% of the CLOUDY solution, except from the warm tail 
extending beyond the location of the ionization front pro- 
duced by CLOUDY94. The nature of this feature is un- 
clear, but it can be probably imputed to heat conduction, 
implemented in CLOUDY94, which cause an heat trans- 
fer across the ionization front where the temperature gra- 
dient is very high. As for the ionization fractions, the errors 
plotted on the CLOUDY94 curves are evaluated as follows. 
At some distances from the source, the CLOUDY94 cal- 
culation gives the unphysical results {xjf+ + Xjfo) > 1 or 
{xue° + ^He+ + ^He++) > Ij ^6 then assumc the errors to 
be the maximum value of the quantities [{xjj+ + Xf^o) — 1], 
for hydrogen species, and [{xjj^o + Xff^+ + Xu^++) — 1] for 
helium species. The errors on the CRASH results are shown by 
the scatter of the plotted points, which is more significant 
for the neutral fractions as a result of the discretization error 
connected to continuous process. The x^o and XHf.o profiles 
produced by the two codes are characterized by the same 
trend. The CRASH results tend to overpredict these quanti- 
ties with respect to CLOUDY94 on average by 20-30%; this 



is also the case for the Xfj^+ , and is likely due to differ- 
ences in the expressions for the rate coefficients adopted. 
For example, at T = 10'' K the CRASH total hydrogen re- 
combination rate (see Appendix B) is equal to 5.097x10"'^ 
cm^ s~\ while that adopted by CLOUDY (see Table 20 
of CLOUDY90 manual) is 4.17x10"^^ cm^ s"', i.e. a 20% 
discrepancy. Our higher value is consistent with the larger 
abundance of H'^ found inside the Hn. The Xfi+ profiles are 
in perfect agreement, as well as the X[j^++ profiles which 
appear to be in very good agreement within the errors. 

To get more insights on the numerical technique 
adopted, we plot in Fig. |S] the evolution of the temperature 
and ionization fractions for the various species in four fixed 
cells at different distances from the source {d = 4, 14, 40 and 
52 pc). It is important to note that the continuity of the 
photoionization (photoheating) process is well reproduced 
at all distances; also, cells at different distances start to be 
ionized at different times, as the I-front propagates into the 
density field. Only in the cell closest to the source the radia- 
tion field at energies higher than 4 Ryd is strong enough to 
significantly doubly ionize He; at 14 pc from the source He"'' 
is only partially ionized, and in the two most distant cells 
He"*"^ aboundance is almost negligible. The flux depletion 
and the progressive filtering of the spectrum by the inter- 
vening medium, are further illustrated by Fig. |S| where we 
plot the frequency distribution of photon packets crossing 
the same four cells as above during the entire simulation, 
weighted by their photon content. Moving to regions farther 
from the source, the ionizing flux decreases, as expected, due 
to the combined effect of geometrical dilution and photo- 
electric absorption; this last process also causes the filtering 
effect resulting in the progressive hardening of the spectrum. 
Again, we note that packets with energy above 4 Ryd are 
heavily depleted in distant regions. 

These tests allow us to conclude that CRASH reliably 
computes, on limited computational times (the discussed 
simulation took a few hours on a 1 GHz workstation), the 
physical quantities of a photoionized gas of primordial com- 
position and, when confronted to state-of-the-art 1-D equi- 
librium codes as CLOUDY94, produces results of compara- 
ble accuracy level. 

As a final remark, we note that the presented tests con- 
stitute a particularly difficult physical benchmark due to the 
relatively high densities and soft source spectrum. 

3.3 Multiple Point Sources 

As a third test we run the case of multiple radiation sources 
in the box. Again, we consider a homogeneous H/He den- 
sity field with n = 1 cm"^, Nc = 128, Np = 10* and 
Lbox ~ 160 pc. Two sources are embedded in the compu- 
tational volume: one with Li = 10^* erg s~^, and a weaker 
one with L2 = Li/27, both with a T = 60000 K black-body 
spectrum. The temperature and ionization fractions fields 
are initialized as in the previous test. We follow the evolu- 
tion of this system for a time ts = 4.5 x 10^ yr « S-it^ec- 

In Fig. 13 we show the ionized hydrogen fraction distri- 
bution in a slice of the simulation box through the plane 
containing the source locations. The panels refer to four dif- 
ferent simulation times: (0.3,0.6,0.9,3) x t^f.^. The Hn re- 
gions around each source evolve independently until they 
merge at t ~ t^ec- Once the Hn regions overlap, photons 



10 A. Maselli, A.Ferrara and B.Ciardi 




from one source can change the properties of the ionized re- 
gion previously produced by the other one. In particular, we 
notice that the I-front curvature of the smaller Hn region 
is changed in the overlapping region. This effect is better 
illustrated by Fig. |S] (left panel), where a 1-D cut of the 
ionized hydrogen fraction along a direction parallel (« 10 
cells away) to that connecting two sources is shown. The 
three lines refer to simulations in which either one or the 
other source were turned off and in which both sources are 
active. As described above, in the cumulative curve an en- 
hancement of the fraction is clearly visible both in the 
overlapping region and on the far side of the smaller Hn re- 
gion due to percolating photons coming from the stronger 
source. A similar effect is visible in the temperature profiles 
along the line of sight connecting the two sources, plotted 
in the right panel of Fig. |H| In the overlapping region the 
temperature is larger than for the case of a single source. 

In order to understand the reason for the difference in 
the temperature profiles within the Hn regions of the two 
sources we have performed an additional check. The heating 



rate associated with photoionizations of the species A can be 
written as T-La ~ VA^hv), where (Jiv) is the mean energy of 
the photoelectrons. As the spectrum is the same, TIa only 
depends on the photoionization rate. The temperature re- 
sults from the balance of such photoheating with radiative 
losses. In turn, we have checked that the latter are largely 
dominated by the recombination cooling rate, which can also 
be written as the product of the recombination rate and the 
mean energy lost per recombining electron. Fig.|^shows the 
time evolution of the temperature, hydrogen ionization frac- 
tion (upper panel) , photoionization and recombination rates 
(lower panel) for two cells located at a distance d = 5 pc 
from each source. The different luminosity of the sources re- 
sults in a different gas thermal/ionization history, the cell 
illuminated by the most luminous source reaching equilib- 
rium earlier. The photoionization rate increases smoothly 
with time, due to the decrease of the intervening opacity, 
until it reaches a maximum; thereafter it decreases because 
of the lack of absorbers. The temperature equilibrium value 
is reached in both cases slightly before the recombination 
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Figure 5. Evolution of the temperature, T (in units of 10* K; solid lines) and ionization fractions for x-^j^ (dashed), (dotted-dashed) 
and Xu^-{-+ (long dashed), in four fixed cells at different distances from the source, d=l, 14, 40 and 52 pc. 
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Figure 6. Frequency distribution of the photon packets, weighted by their photon content, seen by the same four cells of Fig.l^over the 
entire simulation. 
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Figure 7. Slice of tlie simulation box through the planes containing the source locations showing the ionized hydrogen fraction, , 
distribution at four different simulation times: (0.3, 0.6, 0.9, 3) X t!j!^^ (from left to right and top to bottom). 
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Figure 8. Left: Ionized hydrogen fraction along a direction parallel (fs; 10 cells away) to the one connecting the two sources and passing 
trough the I-fronts intersection region. The three lines refer to simulations in which either the stronger (dashed line) or the weaker source 
(dotted-dashed) were turned off; the third cut (solid) is through the simulations in which both sources arc active. Right: Temperature 
along the line connecting the two sources; the notation for the curves is the same as in the left panel. 
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Figure 9. Upper panel: Time evolution of the temperature (solid), hydrogen ionization fraction (dashed) for two cells located at a 
distance d = 5 pc from each source. The label 1 (2) on the curves refer to the cell at 5 pc from the the source wiht luminosity Li 
{L2=Li/27). Lower panel: As above, but for the photoionization (solid) and recombination (dashed) rates. 



rate reaches a plateau at « 10~^ yr~^; this value is the 
same for the two sources as it essentially depends only on 
the gas density (the tiny deviation seen at large evolutionary 
times is due to the temperature dependence of the recom- 
bination coefficient). However, the photoionization rate at 
that stage is larger (by roughly a factor Li/L2=27) for the 
more luminous source than for the fainter one. This extra 
net heating explains the temperature difference. It is also 
worth noting that asymptotically the photoionization and 
the recombination rates become equal as expected from the 
equilibrium condition for a Stromgren sphere. 

Finally, in Fig.[Tniwe plot 2-D maps of the , x^^+ , 
^He++ temperature distributions at the end of the sim- 
ulation. Due to the larger mean free path typical of He" 
ionizing photons, the He"^ I-front is wider than the H"*" one, 
and it is also smoother due the flatter spectrum at hi/ ^ 24.6 
eV. Moreover, the He""" I-front appears to be more affected 
by numerical noise because of the poorer sampling of the 
helium ionizing region of the spectrum. This is even more 
evident in the He"""^ distribution for which the I-front is 
jagged in the outer parts, where ionization balance is gov- 
erned by the very high energy tail of the power spectrum. 
Nevertheless the calculation correctly reproduces the higher 
He"*" ionization level around to the stronger source. 



3.4 Shadowing of background radiation 

As a final test, we study the case of a gas ionized by a 
power-law background radiation field. Our main aim is to 
demonstrate the ability of CRASH to deal with the problem of 
diffuse recombination radiation. The density field has been 
initialized with a high density region (n=l cm~^) of dimen- 
sion 128x100x20 cells, embedded in an homogeneous gas 
with n — 10~^ cm ~^ . The density field is exposed to a 
plane parallel ionization front traveling orthogonally to the 
overdense slice. The background spectrum has a power-law 
shape, J„ = Joiv /vth,Ho)~°' , with Jo = 10~^^ erg cm"^ s"^ 
Hz~^ and a = —1.4. The gas composition is the same as in 
the previous tests. The gas is initially completely neutral and 
at a temperature T = 10^ K in the entire simulation box; 
the physical simulation time is = 3 x 10* yr« A.'it^^^, the 
box linear size is Lbox ~ 6.6 kpc; moreover, Np = 5 x lO'^ 
and Nc = 128. 

Fig. 1101 shows maps of the ionized hydrogen fraction distri- 
bution, orthogonal to the overdense slice, at four different 
simulation times: (0.2, 0.4, 0.9, 2) x t^?^^. The I-front initially 
propagates into the low density field, but it is stopped at the 
edge of the slice by the high recombination rate, producing a 
shadow behind it. The I-front is very smooth due the power- 
law spectrum of the ionizing background. As the simulation 
proceeds, the diffuse radiation produced by recombinations 
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Figure 10. Slices of the simulation box through the plane containing the sources locations showing Xjj+ , , Xjj^++ and temperature 

distributions (from left to right and from top to bottom) at the end of the simulation, t = 3At^^^. 



in the ionized gas, begin to ionize the low density region be- 
hind the high density shce. The diffuse radiation then prop- 
agates into the shadowed region, ionizing gas that would 
otherwise remain neutral. Finally, Fig. 1111 shows maps for 
the final X[j+, Xu^+, X[j^++ and temperature distributions. 
The hard spectrum is able to ionize almost completely he- 
lium to He"^"'" in the regions exposed to direct radiation, 
whereas the diffuse radiation is not hard enough to produce 
the same effect. In the He^"'" distribution (lower- left panel), 
a clear separation between the gas ionized by either direct 
or diffuse radiation is visible. A similar signature is present 
also in the temperature map, again caused by the fact that 
diffuse photons have typical energies confined in a narrow 
band above ionization thresholds of the three absorbers; as 
a consequence, little energy is available to photoheat the 
gas. In conclusion, although necessarily qualitative due to 
the lack of precise solutions to this problem, the shadowing 
experiment performed here seems to yield physically plausi- 
ble results. 



4 SUMMARY 

We have presented an updated version of CRASH , a 3-D RT 
code evaluating the effects of an ionizing radiation field prop- 
agating through a given inhomogeneous H/He density field 
on the physical conditions of the gas. The code, based on a 



Monte Carlo technique, self-consistently calculates the time 
evolution of gas temperature and ionization fractions due to 
an arbitrary number of point/extended sources and/or dif- 
fuse background radiation with given spectra. In addition, 
the effects of diffuse ionizing radiation following recombi- 
nations of ionized atoms have been included. The code has 
been primarily developed to study a number of cosmological 
problems, such as hydrogen and helium reionization, phys- 
ical state of the Lya forest, escape fraction of Lyman con- 
tinuum photons from galaxies, diffuse Lya emission from 
recombining gas, etc. However, its fiexibility allows appli- 
cations that could be relevant to a wide range of astro- 
physical problems. The code architecture is sufficiently sim- 
ple that additional physics can be easily added using the 
same algorithms already described in this paper. For exam- 
ple, dust absorption/re-emission can be included with min- 
imum effort; molecular opacity and line emission, although 
more complicated, do not represent a particular challenge 
given the numerical scheme adopted. Obviously, were such 
processes added, the computational time could become so 
long that parallelization would be necessary. This would 
be required also when CRASH will be coupled to a hydro- 
dynamical code to study the feedback of photo-processes 
onto the (thermo-)dynamics of the system. This perspective 
development looks quite encouraging: in fact, Monte Carlo 
schemes are by construction extremely suitable for paral- 
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lei computations, as the packet load can be distributed in 
a straightforward manner to processors, limiting the inter- 
communication problems among them to a very minimum. 
Clearly, this is a noticeable advantage of the adopted scheme 
over those based on full solution of the RT equations. This 
study is already under work and its implementation will be 
presented in a forthcoming communication. 

To demonstrate the performances, accuracy and 
robustness of the code we have studied in detail four 
different test cases designed to ascertain specific aspects 
of RT: (i) pure hydrogen isothermal Stromgren spheres; 
(ii) realistic Stromgren spheres; (iii) multiple overlapping 
point sources, and (iv) shadowing of background radiation. 
When possible, a detailed quantitative comparison of the 
results against either analytical solutions or 1-D standard 
photoionization codes has been made showing a remarkable 
level of agreement ( ^ few percent). For more challenging 
tests the code yields physically plausible results, which 
could be checked only by comparison with other similar 
codes. 

This work was partially supported by the Research and 
Training Network 'The Physics of the Intergalactic Medium' 
set up by the European Community under the contract 
HPRN-CT2000-00126 RG29185. 
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APPENDIX A: MONTE CARLO SAMPLING 

The Monte Carlo (MC) method, extensively used in many 
field of scientific research, allows to sample with great effi- 
ciency physical quantities with a given Probability Distribu- 
tion Function (PDF). Here we outline the basic principles on 



which the MC sampling technique is based. Let us consider 
a variable q having values in the interval [a,b], representing 
a physical quantity distributed according to a PDF, /(g), 
normalized to J^^ f{<l)dq — 1. The PDF is defined so that 
the probability of having a value q in the range [q' , q' + dq] 
is f{q')dq. Let us also introduce the cumulative probability 
P{q'), defined as: 



f{<l)dq, 



(Al) 



which measures the probability that the previous value of 
q will yield a value lower or equal to q' . Given the above 
normalization, P{q') ^ 1 for each q' in [a,b]. 

Once the PDF is given, a number TZ is randomly ex- 
tracted in the interval [0, 1]. Then, q' is derived setting 
TZ — P{q') and inverting the integral in eg. lAll In this way 
the probability of obtaining the value q' when extracting 
randomly a number in [0,1], corresponds exactly to f{q')- 

If the number of extractions is large enough, it is possi- 
ble to obtain a statistical population of values which repro- 
duces with great accuracy the given PDF. The main advan- 
tage of this technique lies in the ability to sample arbitrary 
distribution functions by simply randomly sampling the in- 
terval [0, 1], an easy and cheap algorithm to implement. 
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He° : VHeo{T) = 1.55xlO"^®T' 



26Tn0.3647 



APPENDIX B: RATE COEFFICIENTS AND • Recombination cooling [erg cm^ s"^] (2) 

CROSS SECTIONS 

In this Section we list the rate coefficients and the cross ^rp\ — g yoxlO^^^T^^^ ( ^ \ 

sections adopted by CRASH , for all the physical processes ' ^ ' ' \ iqs J 

included. 

• Photoionization cross sections [cm~^] (1) 

H° : a«oM = 6.3xlO-i**(/.M^,Ho)"' 
He° : cT«,o(z/) = 7.2xlO-^«[1.66(i./i.,ft,Heo)"' °' 

+0.66(z//!/,;,,Hao)-'°'] 

He+ : OHe+iv) = 1.58 x 10-'*(i./z/t^,He+)"^ 

• Recombination rates [cm^ s~^] (2) 



He+ : VHe+{T) = 3.48x IQ-^^T^/^ (3^) ° ' 



Collisional excitation cooling (2) 



He° : aH,o{T) = 1.50x 10-'°r-°-'^'^ 
He+ : ane+{T) = 3.36xlO"^°T 



H° : Vifo(r) = 7.5x10"^" 

[V'/fo]=[erg cm^ s"^] 
He+ : i^He+iT) = 5.54x10-''^ T-" '^^'^ 
-1 ['/'ife+l=[ei'g cm^ s-i] 

He° : Vffe(T) = 9.10x 10"^''T-°-^^*'^ 
[tpHe] = [erg cm'' s^^] 



y \i/2' 



-118348/r 



^+(10^) 



y N 1/2 



-473638/T 



-13179/T 



/ T \ -0-2 
lOy-1/2 / _L \ 



• Collisional ionization rates [cm^ s ^] (2) 
H° : 7jfo(T) = 5.85x10""T^/^ 

He° : 7jyeo(r) = 2.38x 10-"T^/^ 



^+(10^) 



He+ : 7Hc+(r) = 5.68x IQ-^^T^/^ 



1 + 



r ^^1/2' 



- 157809. l/T 



-285335.4/T 



-631515/T 



The cooling function, A(T, Ue , n^o , , n^^o , n^(,+ , n^g++ ) , 
is calculated including the contribution of the following 1981; (4) Haiman et al. 1996 
radiative processes: 

• Collisional ionization cooling [erg cm^ s~^] (2) 



The above expressions are for excitations to all levels, to 
n = 2 for He+ and to the n = 2, 3, 4 He° triplets (of the He°(23S) 
state, supposed to be populated by He"*" recombinations) 

• Bremsstrahlung cooling [erg cm~^ s~^] (3) 
/3(T) = 1.42 X lO-'^r^/^ j^^^ +nHe+ +4n«,++lne 

• Compton cooling/heating [erg cm~^ s~^] (4) 
w{T) = 1.017 X lO'^'^T* [T - T^] rie 

References: (1) Osterbrok 1989; (2) Cen 1992; (3) Black 



H° : Cifo(r) = 1.27x10"^ V^/^ 
He° : Cifeo(r) = 9.38xl0"^^r^/^ 
He+ : Cne+iT) = 4.95x 10"^^T^/^ 



rp .1/2- 



1 + 



rp . 1/2' 



10^ 



J, .1/2' 



VlOV 



-157809. l/T 



-285335.4/T 



-631515/T 



